Computer-aided modeling and manufacturing requires precise control of spatial position along a prescribed motion path. We will illustrate the use of Adaptive Quadrature to solve a fundamental piece of the problem: equipartition, or the division of an arbitrary path into equal-length subpaths.
In numerical machining problems, it is preferable to maintain constant speed along the path. During each second, progress should be made along an equal length of the machine– material interface. In other motion planning applications, including computer animation, more complicated progress curves may be required: A hand reaching for a doorknob might begin and end with low velocity and have higher velocity in between. Robotics and virtual reality applications require the construction of parametrized curves and surfaces to be navigated. Building a table of small equal increments in path distance is often a necessary first step.
Assume that a parametric path P = $\{x(t),\ y(t) \ | \ 0 \leq t \leq 1\}$ is given. Figure 5.6 shows the example path
$$ P = \begin{cases} x(t) = 0.5 + 0.3t + 3.9t^2 − 4.7t^3, \\ y(t) = 1.5 + 0.3t + 0.9t^2 − 2.7t^3 \end{cases} $$
Bézier curve defined by the four points $(0.5,1.5)$, $(0.6,1.6)$, $(2,2)$, $(0,0)$. Points defined by evenly spaced parameter values $t = 0, 1/4, 1/2, 3/4, 1$ are shown.
Note that even spacing in parameter does not imply even spacing in arc length. Your goal is to apply quadrature methods to divide this path into $n$ equal lengths.
Recall from calculus that the arc length of the path from $t_1$ to $t_2$ is
$$ \int_{t_1}^{t_2} \sqrt{x'(t)^2+y'(t)^2}dt $$
Only rarely does the integral yield a closed-form expression, and normally an Adaptive Quadrature technique is used to control the parametrization of the path.
Commence by defining global imports, the given parametric functions and its curve, as well as a global tolerance and the number of decimal places that will be specified throughout the project.
""" Global Imports """
# Plots and animation
import matplotlib
from matplotlib.backends.backend_agg import FigureCanvasAgg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation, rc
from celluloid import Camera
import matplotlib.pyplot as plt
from functionsV2 import * # Project Functions. See Appendix.
from time import perf_counter # For timing
import numpy as np
from random import uniform
tol = 0.5e-5 # Temporary global tol
dec = 5 # Temporary global number of decimals
# The given parametric equation
x = lambda t: 0.5 + 0.3*t + 3.9*t**2 - 4.7*t**3
y = lambda t: 1.5 + 0.3*t + 0.9*t**2 - 2.7*t**3
""" Plot of the given parametric path """
plotParameterizedCurves(x,y,[0,1],1,"Original Parametrized Curve")
Write a Python function that uses Adaptive Quadrature to compute the arc length from $t_1 = 0$ to $t_2 = T$ for a given $T \leq 1$.
""" Suggested Activity 1 """
# Manually computed derivatives of the given parametric path
dxdt = lambda t: 0.3 + 7.8*t - 14.1*t**2
dydt = lambda t: 0.3 + 1.8*t - 8.1*t**2
# Compute the lambda function that is to be integrated
f = sqrtFunSquared(dxdt, dydt)
# Compute the corresponding arc length
# by computing the integral of f from
# 0 to 1 using the method of Adaptive Quadrature
arcLength = compArcLength(f, 0.0, 1.0, adQuad, tol)
print("The arc length is", round(arcLength, dec))
Write a program that, for any input $s$ between $0$ and $1$, finds the parameter $t^∗(s)$ that is $s$ of the way along the curve. In other words, the arc length from $t = 0$ to $t = t^∗(s)$ divided by the arc length from $t = 0$ to $t = 1$ should be equal to $s$.
In other other words:
$$ s = \frac{\int_{0}^{t^*(s)} \sqrt{x'(t)^2+y'(t)^2}dt}{ \int_{0}^{1} \sqrt{x'(t)^2+y'(t)^2}dt}$$
Use the Bisection Method to locate the point $t^∗(s)$ to three correct decimal places.
""" Suggested Activity 2 """
# Let's time the function tStarOfS for s = 0.5
s = 0.5
start = perf_counter()
tStar2 = tStarOfSBisect(f, s, adQuad, tol)
end = perf_counter()
elapsedTime2 = end - start
print('The optimal value of t for s =', round(s, dec), 'is', round(tStar2, dec))
print('and was computed in', round(elapsedTime2, dec), 'seconds')
# Let's verify if tStar really is the root of the function
print('Which is the root within our tolerance:',
(np.abs(compArcLength(f, 0.0, tStar2, adQuad, tol) /
arcLength - s) < 2*tol))
Equipartition (split into parts of equal length) the path of Figure 5.6 into $n$ subpaths of equal length, for $n = 4$ and $n = 20$. Plot analogues of Figure 5.6, showing the equipartitions. If your computations are too slow, consider speeding up the Adaptive Quadrature with Simpson’s Rule, as suggested in Computer Problem 5.4.2.
Let us consider the scenario where $n = 4$:
""" Suggested Activity 3 """
# For n = 4
sArray_n4 = [0.0, 0.25, 0.5, 0.75, 1.0]
n = 4
tStarArray = np.zeros(n+1)
start = perf_counter()
for i in range(n):
tStarArray[i+1] = tStarOfSBisect(f, sArray_n4[i+1], adQuadSimpson, tol)
end = perf_counter()
elapsedTime3n4 = end - start
print("For n = 4 the computation of t*(s) took",
round(elapsedTime3n4, dec), "seconds using Simpson's Rule and the Bisection method")
for i in range(n): # let's verify that each arclength is about a quarter of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
# Plot the equipartitioned subpaths on one graph but in different colors
plotParameterizedCurves(x,y,tStarArray,n,"Parametrized Curve Equipartitioned into 4 subpaths")
Let us then consider the scenario where $n = 20$:
# For n = 20
sArray_n20 = np.arange(0.00, 1.05, 0.05)
n = 20
tStarArray = np.zeros(n+1)
start = perf_counter()
for i in range(n):
tStarArray[i+1] = tStarOfSBisect(f, sArray_n20[i+1], adQuadSimpson, tol)
end = perf_counter()
elapsedTime3n20 = end - start
print("For n = 20 the computation of t*(s) program took",
round(elapsedTime3n20, dec), "seconds using Simpson's Rule and the Bisection method")
for i in range(n): # let's verify that each arclength is about one-twentieth of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
# Plot the equipartitioned subpaths on one graph but in different colors
plotParameterizedCurves(x,y,tStarArray,n,"Parametrized Curve Equipartitioned into 20 subpaths")
Replace the Bisection Method in Step 2 with Newton’s Method, and repeat Steps 2 and 3.
The derivative with respect to $t$ of the function that is being set to zero, namely $$ d \bigg( s \int_{0}^{1} \sqrt{x'(t)^2+y'(t)^2}dt - \int_{0}^{t^*(s)} \sqrt{x'(t)^2+y'(t)^2}dt \bigg) dt^{-1} $$
The derivative can be estimated in the necessary points using Python's Autograd. However, a loop of trials revealed that the Three-point centered-difference formula described in Chapter 5.1 gave the same estimates as Autograd to at least five decimal places, so that will be used instead.
As in Suggested Activity 2, it is clear that the root lies in the interval $[0,1]$. When repeating Suggest Activity 2, it makes sense to use the value of $t^*(s)$ that was computed earlier as an initial guess. That makes Newton's method converge much faster than the Bisection method, as shown below. However, this is an unfair comparison. More precision is however possible when the path is being equipartitioned into $n$ parts. When diving the path into $n$ paths of equal length, one can use the initial guess of $0$ for the first $t^*(s)$, and then one can use that $t^*(s)$ as an initial guess for the next $t^*(s)$, and so on.
Apparently that the initial guess for Newton's method has quite an effect on the time it takes to converge. Our tests showed no clear sign of one method generally being faster than the other, unless the initial guess is extraordinarily good. When performing large-scale computations, it might be good to perform a few iterations of the Bisection method first to get a decent estimate of a starting point for Newton's method, and then proceed from there with Newton's method.
""" Suggested Activity 4 """
start = perf_counter() # Newt takes 12 times longer with initial guess of 0.0 than with 0.8
tStar4 = tStarOfSNewton(f, s, adQuad, 0.8, tol)
end = perf_counter()
elapsedTime4 = end - start
print('The computed value of t*(s) for s =', s, end=' ')
print('using Newton\'s method is', round(tStar4, dec+3))
print('The computed value of t*(s) for s =', s, end=' ')
print('using the Bisection method is', round(tStar2, dec+3))
print('Which is the same within our tolerance:', abs(tStar4-tStar2)<tol, "\n")
print('Time required to compute t*(s) using the Newton\'s method is:',
round(elapsedTime4, dec+2), "seconds")
print('Time required to compute t*(s) using the Bisection method is:',
round(elapsedTime2, dec+2), "seconds")
# For n = 4
n = 4
tStarArray = np.zeros(n+1)
start = perf_counter()
for i in range(n):
tStarArray[i+1] = tStarOfSNewton(f, sArray_n4[i+1], adQuadSimpson, tStarArray[i], tol)
end = perf_counter()
elapsedTime4n4 = end - start
print("For n = 4, using Newton's method, the computation of t*(s) took",
round(elapsedTime4n4, dec), "seconds")
print("For n = 4, using the Bisection method, the computation of t*(s) took",
round(elapsedTime3n4, dec), "seconds")
for i in range(n): # let's verify that each arclength is about a quarter of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
# For n = 20
n = 20
tStarArray = np.zeros(n+1)
start = perf_counter()
for i in range(n):
tStarArray[i+1] = tStarOfSNewton(f, sArray_n20[i+1], adQuadSimpson, tStarArray[i], tol)
end = perf_counter()
elapsedTime4n20 = end - start
print("For n = 20, using Newton's method, the computation of t*(s) took",
round(elapsedTime4n20, dec), "seconds")
print("For n = 20, using the Bisection method, the computation of t*(s) took",
round(elapsedTime3n20, dec), "seconds")
for i in range(n): # let's verify that each arclength is about one-twentieth of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
Use animation commands to demonstrate traveling along the path, first at the original parameter $0 \leq t \leq 1$ speed and then at the (constant) speed given by $t^∗(s)$ for $0 \leq s \leq 1$.
animateBezierCurve(x, y, np.arange(0.0,1.05,0.05),
"Traveling at speed with intervals $dt=0.05$ from $t=0$ to $t=1$")
animateBezierCurve(x, y, tStarArray,
"Traveling at constant speed given by t*(s)")
Experiment with equipartitioning a path of your choice. Build a design, initial, etc. of your choice out of Bézier curves, partition it into equal arc length segments, and animate as in Step 5.
""" Suggested Activity 6 """
n_BezierCurves = 3
BezierCurves, DerivBezierCurves = bezierCurves(n_BezierCurves)
n = 20 # Equipartition each BezierCurve into 20 subpaths
tStarArrays = [np.zeros(n+1),np.zeros(n+1),np.zeros(n+1)]
for j in range(n_BezierCurves):
random_f = sqrtFunSquared(DerivBezierCurves[j][0], DerivBezierCurves[j][1])
random_arcLength = compArcLength(random_f, 0.0, 1.0, adQuadSimpson, tol)
print('The arc length of random bezier curve number', j+1, 'is', random_arcLength)
for i in range(n):
tStarArrays[j][i+1] = tStarOfSNewton(random_f, sArray_n20[i+1], adQuadSimpson, tStarArrays[j][i], tol)
# Verify that each arclength is about one-twentieth of the length of the path
partialArclength = compArcLength(random_f, tStarArrays[j][i], tStarArrays[j][i+1], adQuadSimpson, tol)
if abs(partialArclength / random_arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArrays[j][i], dec), end=' ')
print('to', round(tStarArrays[j][i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / random_arcLength), dec))
plotParameterizedCurves(BezierCurves[j][0],BezierCurves[j][1], tStarArrays[j], n,
"Parametrized Curve Equipartitioned into 20 subpaths",
"RandomBezierCurveEquipartitioned"+str(j))
animateBezierCurve(BezierCurves[0][0], BezierCurves[0][1], np.arange(0.0,1.05,0.05),
"Traveling at speed with intervals $dt=0.05$ from $t=0$ to $t=1$")
animateBezierCurve(BezierCurves[0][0], BezierCurves[0][1], tStarArrays[0],
"Traveling at constant speed given by t*(s)")
animateBezierCurve(BezierCurves[1][0], BezierCurves[1][1], np.arange(0.0,1.05,0.05),
"Traveling at speed with intervals $dt=0.05$ from $t=0$ to $t=1$")
animateBezierCurve(BezierCurves[1][0], BezierCurves[1][1], tStarArrays[1],
"Traveling at constant speed given by t*(s)")
animateBezierCurve(BezierCurves[2][0], BezierCurves[2][1], np.arange(0.0,1.05,0.05),
"Traveling at speed with intervals $dt=0.05$ from $t=0$ to $t=1$")
animateBezierCurve(BezierCurves[2][0], BezierCurves[2][1], tStarArrays[2],
"Traveling at constant speed given by t*(s)")
Write a program that traverses the path P according to an arbitrary progress curve $C(s), 0 ≤ s ≤ 1$, with $C(0) = 0$ and $C(1) = 1$. The object is to move along the curve $C$ in such a way that the proportion $C(s)$ of the path’s total arc length is traversed between $0$ and $s$. For example, constant speed along the path would be represented by $C(s) = s$.
Try progress curves $C(s) = s^{1/3}$, $C(s) = s^2$, $C(s) = sin(s\pi/2)$, for example.
# C(s)=s^2
n = 20
sSquaredArray_n20 = [s**2 for s in sArray_n20]
tStarArray = np.zeros(n+1)
for i in range(n):
tStarArray[i+1] = tStarOfSNewton(f, sSquaredArray_n20[i+1], adQuadSimpson, tStarArray[i], tol)
for i in range(n): # let's verify that NONE of the arclengths is about one-twentieth of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if not abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
plotParameterizedCurves(x,y,tStarArray,n,"Parametrized Curve partitioned into 20 subpaths according $C(s)=s^2$")
animateBezierCurve(x, y, tStarArray,
"Traveling at varied speed given by $C(s)=s^2$")
# C(s)=s^(1/3)
n = 20
sPowOneThirdArray_n20 = [s**(1/3) for s in sArray_n20]
tStarArray = np.zeros(n+1)
for i in range(n):
tStarArray[i+1] = tStarOfSNewton(f, sPowOneThirdArray_n20[i+1], adQuadSimpson, tStarArray[i], tol)
for i in range(n): # let's verify that NONE of the arclengths is about one-twentieth of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if not abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
plotParameterizedCurves(x,y,tStarArray,n,"Parametrized Curve partitioned into 20 subpaths according $C(s)=s^{1/3}$")
animateBezierCurve(x, y, tStarArray,
"Traveling at varied speed given by $C(s)=s^{1/3}$")
# C(s) = sin(s\pi/2)
n = 20
sSinArray_n20 = [np.sin(s*np.pi/2) for s in sArray_n20]
tStarArray = np.zeros(n+1)
for i in range(n):
tStarArray[i+1] = tStarOfSNewton(f, sSinArray_n20[i+1], adQuadSimpson, tStarArray[i], tol)
for i in range(n): # let's verify that NONE of the arclengths is about one-twentieth of the length of the path
partialArclength = compArcLength(f, tStarArray[i], tStarArray[i+1], adQuadSimpson, tol)
if not abs(partialArclength / arcLength - 1/n) > 10*tol:
print('Arc length from', round(tStarArray[i], dec), end=' ')
print('to', round(tStarArray[i+1], dec), end=' ')
print('is', round(partialArclength, dec))
print('Proportional arc length:',
round(np.abs(partialArclength / arcLength), dec))
plotParameterizedCurves(x,y,tStarArray,n,"Parametrized Curve partitioned into 20 subpaths according $C(s) = sin(s\pi/2)$")
animateBezierCurve(x, y, tStarArray,
"Traveling at varied speed given by $C(s)=s^{1/3}$")